#include "TimedVoxels.hpp"

#include <iostream>

#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <cmath>
#include <algorithm>

#include "config.hpp"
#include "misc.hpp"

using namespace std;

namespace Dynamic
{

static float image[196 * 128], image1[196 * 128];

int TimedVoxels::VoxelIndex(int i, int j)
{
	return (i * vrows + j);
}

TimedVoxels::TimedVoxels(int timeSliceCount, double dt) : timeSliceCount(timeSliceCount), dt(dt)
{
	this->array = new double[nvoxels * this->timeSliceCount];
}

TimedVoxels::TimedVoxels(const TimedVoxels& v)
{
	this->array = nullptr;
	*this = v;
}

TimedVoxels::~TimedVoxels()
{
	delete [] array;
}

double& TimedVoxels::at(int i, int j, unsigned int timeSliceId)
{
	if (i < 0 || i > vrows - 1 || j < 0 || j > vrows - 1) {
		throw std::exception("Index out of range");
	}
	return this->at(VoxelIndex(i, j), timeSliceId);
}

const double& TimedVoxels::at(int i, int j, unsigned int timeSliceId) const
{
	if (i < 0 || i > vrows - 1 || j < 0 || j > vrows - 1) {
		throw std::exception("Index out of range");
	}
	return this->at(VoxelIndex(i, j), timeSliceId);
}

double& TimedVoxels::at(int voxelIndex, unsigned int timeSliceId)
{
	return array[voxelIndex * timeSliceCount + timeSliceId];
}

const double& TimedVoxels::at(int voxelIndex, unsigned int timeSliceId) const
{
	return array[voxelIndex * timeSliceCount + timeSliceId];
}

void TimedVoxels::operator=(const TimedVoxels& m)
{
	if (&m != this) {
		this->timeSliceCount = m.timeSliceCount;
		this->dt = m.dt;

		delete[] this->array;
		this->array = new double[nvoxels * this->timeSliceCount];
		for (unsigned int i = 0; i < nvoxels * this->timeSliceCount; i++) {
			this->array[i] = m.array[i];
		}
	}
}

void TimedVoxels::Zero()
{
	for (unsigned int V = 0; V < nvoxels * this->timeSliceCount; ++V) {
		array[V] = 0;
	}
}

void TimedVoxels::start(double hits)
{
	// Spread the expected emission density over the voxels (in space and time)
	for (unsigned int V = 0; V < nvoxels * this->timeSliceCount; ++V) {
		array[V] = (double)hits / (nvoxels * this->timeSliceCount);
	}
}

double TimedVoxels::getMeasurementInterval() const
{
	return this->dt * this->timeSliceCount;
}

void TimedVoxels::InitVoxels(double activity /*= 1*/)
{
	double measurementInterval = this->getMeasurementInterval();
	for (unsigned int i = 0; i < vrows; ++i) {
		for (unsigned int j = 0; j < vrows; ++j) {
			for (unsigned int t = 0; t < this->timeSliceCount; ++t) {
				at(i, j, t) = 0;
				if (i >= vrows / 4 && i <= vrows / 4 + 1 && j >= vrows / 4 && j <= vrows / 4 + 1) {
					at(i, j, t) = 16.0;
				}
				if (i >= vrows / 2 && i < vrows / 2 + 7 && j > vrows / 2 && j < vrows / 2 + 7) {
					at(i, j, t) = 1.0;
				}
				at(i, j, t) *= activity * Gauss1D(t * this->dt, measurementInterval * 0.25);
			}
		}
	}
}

double TimedVoxels::getAccumulatedValue(int voxelIndex) const
{
	double vSum = 0;
	for (unsigned int t = 0; t < this->timeSliceCount; ++t) {
		vSum += this->at(voxelIndex, t) * this->dt;
	}

	return vSum;
}

double TimedVoxels::getAccumulatedValue(int i, int j) const
{
	return this->getAccumulatedValue(VoxelIndex(i, j));
}

double TimedVoxels::getValueAtTime(int voxelIndex, double t) const
{
	double sum = 0;
	// Accumulate for all basis functions
	for (unsigned int i = 0; i < this->timeSliceCount; i++) {
		sum += this->at(voxelIndex, i) * linearUniformTentBasisFunction(this->dt, i, t);
	}

	return sum;
}

double TimedVoxels::MaxAccumulated() const
{
	double maxvoxel = 0;
	for (unsigned int V = 0; V < nvoxels; ++V) {
		double vSum = this->getAccumulatedValue(V);
		if (maxvoxel < vSum) {
			maxvoxel = vSum;
		}
	}
	return maxvoxel;
}

double TimedVoxels::Max(int border /*= 0*/) const
{
	double maxvoxel = 0;
	for (unsigned int i = 0 + border; i < vrows - border; ++i) {
		for (unsigned int j = 0 + border; j < vrows - border; ++j) {
			for (unsigned int t = 0; t < this->timeSliceCount; ++t) {
				double val = this->at(i, j, t);
				if (maxvoxel < val) {
					maxvoxel = val;
				}
			}
		}
	}
	return maxvoxel;
}

double TimedVoxels::Sum() const
{
	double sum = 0.0;
	for (unsigned int V = 0; V < nvoxels; ++V) {
		double vSum = this->getAccumulatedValue(V);
		sum += vSum;
	}
	return sum;
}

double TimedVoxels::L2ErrorSum(const TimedVoxels& xRef, double corr /*= 1.0*/) const
{
	double error = 0, sum = 0;
	for (unsigned int i = 0; i < nvoxels; ++i) {
		double vSum = this->getAccumulatedValue(i) * corr;
		double vSumRef = xRef.getAccumulatedValue(i);
		error += sqr(vSumRef - vSum);
		sum += sqr(vSumRef);
	}
	return sqrt(error / sum) * 100;
}

double TimedVoxels::L2ErrorTime(const TimedVoxels& xRef, double corr /*= 1.0*/) const
{
	double error = 0, sum = 0;
	for (unsigned int i = 0; i < nvoxels; ++i) {
		for (unsigned int tId = 0; tId < this->timeSliceCount; ++tId) {
			double t = tId * this->dt;
			double vSum = this->getValueAtTime(i, t) * corr;
			double vSumRef = xRef.getValueAtTime(i, t);
			error += sqr(vSumRef - vSum);
			sum += sqr(vSumRef);
		}
	}
	return sqrt(error / sum) * 100;
}


void TimedVoxels::PrintImage(const char* ofname, double maxvoxel, bool bar /*= false*/)
{
	cout << "\n Image: " << ofname << '\n';
	static char buffer[80];
	sprintf(buffer, "%s.ima", ofname);
	PrintPlot(buffer, maxvoxel);

	sprintf(buffer, "%s.tga", ofname);
	FILE* ofile = fopen(buffer, "wb");
	int Xres = (bar) ? xres : yres;
	PrintTGAHeader(ofile, Xres, yres);

	for (unsigned int i = 0; i < yres; ++i) {
		for (unsigned int j = 0; j < yres; ++j) {
			image[i * vrows * 4 + j] = this->getAccumulatedValue(i / image_scale, j / image_scale) / maxvoxel;
		}
	}
	glClear(GL_COLOR_BUFFER_BIT);
	glRasterPos2f(-1, -1);
	glDrawPixels(yres, yres, GL_LUMINANCE, GL_FLOAT, image);
	sprintf(buffer, "%d", 0);
	Message(buffer, false);
	sprintf(buffer, "%d", (int)maxvoxel);
	Message(buffer, true);
	glRasterPos2f(-1, -1);
	glReadPixels(0, 0, Xres, yres, GL_LUMINANCE, GL_FLOAT, image);


	for (int i =  yres - 1; i >= 0; --i) {
		for (unsigned int j = 0; j < Xres; ++j) {
			if (bar && j >= yres - 3 && j < yres) {
				image[i * Xres + j] = (double)(i - 1) / yres;
			}
			PrintTGAColor(ofile, (double)image[i * Xres + j], 1);
		}
	}
	fclose(ofile);
}

void TimedVoxels::PrintPlot(const char* ofname, double maxvoxel)
{
	cout << "\n Image: " << ofname << '\n';
	FILE* ofile = fopen(ofname, "w");

	for (unsigned int i = 0; i < vrows; ++i) {
		for (unsigned int j = 0; j < vrows; ++j) {
			// Commented out by Laci :)
			//              fprintf(ofile, "%d %d %f\n", i, j, at(i, j) / maxvoxel);
			if (i == vrows / 4 && j == vrows / 4) {
				fprintf(ofile, "%f ", maxvoxel);
			} else {
				double val = this->getAccumulatedValue(i, j);
				fprintf(ofile, "%f ", val < maxvoxel ? val : maxvoxel);
			}
		}
		fprintf(ofile, "\n");
	}
	fclose(ofile);
}

void TimedVoxels::PrintPlotForDensityOverTime(const char* ofname, int minX, int minY, int maxX, int maxY, double dt)
{
	cout << "\n GNUPlot text file (density over time): " << ofname << '\n';
	FILE* ofile = fopen(ofname, "w");

	double measurementTime = this->dt * this->timeSliceCount;

	for (double t = 0.0; t < measurementTime; t += dt) {
		double sum = 0.0;
		for (unsigned int i = 0; i < vrows; ++i) {
			for (unsigned int j = 0; j < vrows; ++j) {
				if (minX <= i && i <= maxX && minY <= j && j <= maxY) {
					sum += getValueAtTime(VoxelIndex(i, j), t);
				}
			}
		}
		fprintf(ofile, "%f %f", t, sum);
		fprintf(ofile, "\n");
	}
	fclose(ofile);
}

void TimedVoxels::PrintInTxt(const char* ofname, double maxvoxel)
{
	cout << "\n Mathematica text file (density over time): " << ofname << '\n';
	FILE* ofile = fopen(ofname, "w");

	int border = 5;
	double measurementTime = this->dt * this->timeSliceCount;
	double maxValue = this->Max(border);
	cout << "maxValue: " << maxValue << '\n';
	cout << "maxvoxel: " << maxvoxel << '\n';

	fprintf(ofile, "{");
	for (unsigned int i = 0 + border; i < vrows - border; ++i) {
		if (i > 0 + border) {
			fprintf(ofile, ", ");
		}
		fprintf(ofile, "{");
		for (unsigned int j = 0 + border; j < vrows - border; ++j) {
			if (j > 0 + border) {
				fprintf(ofile, ", ");
			}
			fprintf(ofile, "{");
			for (unsigned int t = 0; t < this->timeSliceCount; ++t) {
				double val = this->at(i, j, t) / maxValue;
				if (t > 0) {
					fprintf(ofile, ", ");
				}
				fprintf(ofile, "%f", std::min<double>(val, 1.0));
			}
			fprintf(ofile, "}");
		}
		fprintf(ofile, "}");
	}
	fprintf(ofile, "}");
	fclose(ofile);
}


}